% function ret_4=bootstrap_final(data,rep,p,h,var,quan,k);
% p=6;rep=100;
% var= variable for which irf bootstrap is computed;;
% quan=lower pefcentile of confidence interval
% h=num of horizons
% n=num of domestic variables
% k=num of variables;
% m=num of exogenous variables
% rep=num of bootstrap replications;
function [dat,F,c_tot]=boot_us_gen_var_data(data,p,m,trend_us,rest,st,LE,r1,r2,lr,exo,exo_for,quadr,kil,bias_corr_F,bias_corr_C)
% clear;data=randn(200,6);
% m=5;p=5;trend_us=1;

 T=length(data);rep=1 ;h=60;k=length(data(1,:));trend_uk=trend_us;
bayes=0;
%  rest=1;st=4;
  ret=nad(data,0,m,T,p,trend_us,r1,LE,exo,rest,st,bayes,quadr);
sig_x=ret.cov_res;
 coint_x=ret.coint;
 c1=ret.c;
 E=ret.G;
res_x=ret.res;
if LE>0&&trend_us>0
t_x=ret.trend_coef;
end
if quadr==1
t_xx=ret.trend_coef1;
end
SIG1=ret.SIG;
% THE ZERO BELOW RIGHT TO LE IS EXO WHICH IS ALMOST ALWAYS ZERO UNLES YOU WANT FOREIGN VARIABLES TO BE EXOGENOUS WITH NO LAGS%
ret=nad(data,m,k,T,p,trend_uk,r2,LE,exo_for,rest,st,0,quadr);
G=ret.G;
 coint_y=ret.coint;
sig_y=ret.cov_res;
c2=ret.c;
res_y=ret.res;SIG2=ret.SIG;
if trend_uk==1
t_y=ret.trend_coef;
end
if quadr==1
t_yy=ret.trend_coef1;
end
 ret=irf_nad(c1,c2,E,G,r1,r2,sig_x,sig_y,res_x,res_y,coint_x,coint_y,trend_uk,trend_us,m,k-m,T,p,h,LE,lr,rest,st,SIG1,SIG2,bayes,quadr,0,[],[]);
  res_red_y=ret.res_red_y;
  if bayes==1
 ret2.SIGMA=[[ret.SIG;zeros(1,k)] zeros(k,1)];
 end
 
% s_us=ret.shocks_us;
% s_open=ret.shocks;
% TR=repmat(tr',T-p,1);
% TR_=[TR zeros(T-p,(p-1)*k)];
% if LE>0
%     p=p;
% else
%     p=p;
% end
zero2=zeros(k);
% if trend_us<1
%     ret.t_x=zero2(1:m,1);
% else
%     ret.t_x=ret.t_x;
% end
% if trend_uk<1
%     ret.t_red_y=zero2(1:k-m,1);
% else
%     ret.t_red_y=ret.t_red_y;
% end
if LE==1
       R=T-p+2;
   else
       R=T-p+1;
end
   if LE>0&&trend_us>0
tr=[t_x;t_y];
   else
   tr=[zeros(k,1)];
   end
    if quadr==1
trr=[t_xx;t_yy];
    else
         trr=[zeros(k,1)];
    end
TREND=(1:R-1);
TREND_=repmat(TREND,k,1);
TIME=diag(tr)*TREND_;
TREND_FINAL=[TIME;zeros((p-1)*k,R-1)];
TREND1=(1:R-1).^2;
TREND1_=repmat(TREND1,k,1);
TIME1=diag(trr)*TREND_;
TREND_FINAL1=[TIME1;zeros((p-1)*k,R-1)];
if kil==0;
    F=ret.FF;c_tot=ret.c_tot;
else
    F=bias_corr_F;
    c_tot=bias_corr_C;
end
% c_us=ret.IRF_US;
% c_open=ret.IRF_UK;
% d=diag(ret.p);
%  w=s_open*diag(d);
opt = statset('UseParallel','always');
w=[res_x' res_red_y'];
% matlabpool open 4
% [~,u]=bootstrp(rep,@mean,w,'Options',opt);
[~,u]=bootstrp(rep,@mean,w);


% below command is for bootstrapping
res=cell(1,rep);
for i=1:rep
res{i}=[w([u(:,i)],:)];
end
% pack res;
clear u
 for i=2:T
q(:,:,i)=sum(data(1:i,:));
    end
s=squeeze(q(:,:,2:end));
 s=reshape(s',length(data)-1,length(data(1,:)));
 s=[zeros(1,length(data(1,:)));data(1,:);s];
 if LE==1
     data=s(1:end,:);
 end
 if rest==1
             data=[data(:,1:end-st-1) s(2:end,end-st:end)];
 end
for i=1:p
Y{i}=[data(1+p-i:end-i,:)];
end
 Y=cell2mat(Y);
 U=cell(1,rep);
 for i=1:rep
     if LE>0
U{i}=[res{i} zeros(R-1,(p-1)*k)];
     else U{i}=[res{i} zeros(R-1,(p-1)*k)];
     end
 end
 clear res
%  if LE>0
%     p=p+1;;
% else
%     p=p;
% end
    
  C=repmat(c_tot',R-1,1);
  C_=[C zeros(R-1,(p-1)*k)];
  check_data=[F*Y'+[w zeros(R-1,(p-1)*k)]'+C_'+TREND_FINAL+TREND_FINAL1]';
  forecast=[F*Y'+[zeros(R-1,k) zeros(R-1,(p-1)*k)]'+C_'+TREND_FINAL+TREND_FINAL1]';
   test=check_data(:,1:k)-data(p+1:end,:);
   z=cell(1,rep);
for j=1:rep
for i=1:R
if i<=1
    x{i}=Y(i,:)';
else
    x{i}=F*x{i-1}+U{j}(i-1,:)'+C_(i-1,:)'+TREND_FINAL(:,i-1)+TREND_FINAL1(:,i-1);
end
z{j}=x;
end
end
clear U
e=cell(1,rep);
e_=cell(1,rep);
H=cell(1,rep);
for j=1:rep
    e{j}=cell2mat(z{j});
 e_{j}=reshape(e{j}(1:k,:)',R,k);
 H{j}=[data(1:p-1,:);e_{j}];
end
clear z e e_
% H{i}=bootstrapped data
% BELOW COMMANDS IS NORMAL DRAW
% for i=1:rep
% t{i}=randn(T-p,8);
% end
% for i=1:rep
% res{i}=t{i}*diag(d);
% end
% FROM HERE IT'S AS USUAL
% for i=1:rep
%     data_{i}=res{i}-w+data(p+1:end,:);
% end
% for i=1:rep
%     data_{i}=res{i}+data(p+1:end,:);
% end
dat=cell(1,rep);
for i=1:rep
    dat{i}=H{i};
    if LE==1
    dat{i}=[diff(dat{i})];
    else 
        dat{i}=dat{i};
%     data(:,:,i)=dat(:,:,i);
    end
end
    dat=cell2mat(dat);
  if rest>0
        dat=[dat(:,1:k-st-1) [dat(1,k-st:end);diff(dat(1:end,k-st:end))]];
  end
%  dat;clear H
% ret_1=cell(1,rep);sig_x1=cell(1,rep);coint_x1=cell(1,rep);
% E1=cell(1,rep);c11=cell(1,rep);res_x1=cell(1,rep);ret_2=cell(1,rep);ret_3=cell(1,rep);
% G1=cell(1,rep);coint_y1=cell(1,rep);c22=cell(1,rep);sig_y1=cell(1,rep);
% res_y1=cell(1,rep);
%     for i=1:rep
%  ret_1{i}=nad(dat{i},0,m,T,p,trend_us,r1,LE);
% sig_x1{i}=ret_1{i}.cov_res;
%  coint_x1{i}=ret_1{i}.coint;
%  E1{i}=ret_1{i}.G;
%  c11{i}=ret_1{i}.c;
% res_x1{i}=ret_1{i}.res;
%  ret_2{i}=nad(dat{i},m,k,T,p,trend_uk,r2,LE);
% G1{i}=ret_2{i}.G;
%  coint_y1{i}=ret_2{i}.coint;
%  c22{i}=ret_2{i}.c;
% sig_y1{i}=ret_2{i}.cov_res;
% res_y1{i}=ret_2{i}.res;
%  ret_3{i}=irf(c11{i},c22{i},E1{i},G1{i},r1,r2,sig_x1{i},sig_y1{i},res_x1{i},res_y1{i},coint_x1{i},coint_y1{i},trend_uk,trend_us,m,k-m,T,p,h,LE);
%  ret_3{i}.g=[];
%  ret_3{i}.g_=[];
%  ret_3{i}.e=[];
%  ret_3{i}.c_tot=[];
%  ret_3{i}.res_red_y=[];
%  ret_3{i}.SIGMA=[];
%  ret_3{j}.t_x=[];
%  ret_3{i}.t_red_y=[];
%  ret_3{i}.B_1=[];
%  ret_3{i}.B_2=[];
%  ret_3{i}.COINT=[];
%  ret_3{i}.coint_reduced_y=[];
%  ret_3{i}.extended_coint_x=[];
%  ret_3{i}.FF=[];
%  ret_3{i}.IRF=[];
%  ret_3{i}.IRF_US=[];
%  ret_3{i}.IRF_UK=[];
%  ret_3{i}.D=[];
%  ret_3{i}.shocks=[];
%  ret_3{i}.shocks_us=[];
%  ret_3{i}.p=[];
%  ret_3{i}.non_nor_shocks=[];
%   ret_1{i}=[];
%   ret_2{i}=[];
%     end
%      ret_3; ret_1;pack ret_2;pack E1;pack G1;pack coint_x1;
%     pack coint_y1;pack sig_x1;pack sig_y1;pack res_x1;pack res_y1;
%     pack c11;pack c22;pack dat;pack data_;clear ret_1 ret_2 E* G* coint* sig* res* c* s* dat data_
% %     for j=1:rep
% %  ret_3{j}.g=[];
% %  ret_3{j}.g_=[];
% %  ret_3{j}.e=[];
% %  ret_3{j}.c_tot=[];
% %  ret_3{j}.res_red_y=[];
% %  ret_3{j}.SIGMA=[];
% %  ret_3{j}.t_x=[];
% %  ret_3{j}.t_red_y=[];
% %  ret_3{j}.B_1=[];
% %  ret_3{j}.B_2=[];
% %  ret_3{j}.COINT=[];
% %  ret_3{j}.coint_reduced_y=[];
% %  ret_3{j}.extended_coint_x=[];
% %  ret_3{j}.FF=[];
% %  ret_3{j}.IRF=[];
% %  ret_3{j}.IRF_US=[];
% %  ret_3{j}.IRF_UK=[];
% %  ret_3{j}.D=[];
% %  ret_3{j}.shocks=[];
% %  ret_3{j}.shocks_us=[];
% %  ret_3{j}.p=[];
% %  ret_3{j}.non_nor_shocks=[];
% %     end
% 
% %     for j=1:rep
% % RES{j}=[ret_3{j}.S_IRF{1}(1,:)];
% %    end
% % R=reshape(cell2mat(RES'),rep,10);
% % USE COMMANDS BELOW TO ATTAIN THE BOOTSTRAP FOR THE VARIABLE OF INTEREST
% var=5;quan=0.05;
% RES=cell(1,rep);
% for i=1:h+1
% for j=1:rep
% RES{j}=[ret_3{j}.S_IRF{i}(var,:)];
%    end
% R_{i}=reshape(cell2mat(RES'),rep,k);
% end
% pack RES;clear RES
% for i=1:h+1
% Q{i}=[quantile(R_{i},[quan,1-quan])];
% end
% for i=1:h+1
% upper_Q{i}=[Q{i}(2,:)];
% end
% upper_bound=reshape(cell2mat(upper_Q'),h+1,k);
% for i=1:h+1
% lower_Q{i}=[Q{i}(1,:)];
% end
% lower_bound=reshape(cell2mat(lower_Q'),h+1,k);
% % plot([upper_bound(:,5) ret.IRF_US{var}(:,5) lower_bound(:,5)])
% % HALL
% for j=1:k
% for i=1:h+1
% IR{i}=[ret.S_IRF{i}(j:j,1:k)];
% end
% IRF{j}=reshape(cell2mat(IR'),h+1,k);
% end
% lower_bound_h=IRF{var}-(upper_bound-IRF{var});
% upper_bound_h=IRF{var}-(lower_bound-IRF{var});
% % plot([upper_bound_h(:,3) IRF{var}(:,3) lower_bound_h(:,3)])
% % plot([upper_bound_h(:,4) IRF{var}(:,4) lower_bound_h(:,4)])
% % plot([upper_bound_h(:,5) IRF{var}(:,5) lower_bound_h(:,5)])
% % plot([upper_bound_h(:,6) IRF{var}(:,6) lower_bound_h(:,6)])
% plot([upper_bound_h(:,7) IRF{var}(:,7) lower_bound_h(:,7)])
% clear Y TIME TR TREND TREND_ TREND_FINAL IR  Q F C_ C R_ lower_Q lower_bound  r1 r2  tr trend_uk trend_us upper_Q upper_bound  v w x zero2
% ret.g=[];
%  ret.g_=[];
%  ret.e=[];
%  ret.c_tot=[];
%  ret.res_red_y=[];
%  ret.SIGMA=[];
%  ret.t_x=[];
%  ret.t_red_y=[];
%  ret.B_1=[];
%  ret.B_2=[];
%  ret.COINT=[];
%  ret.coint_reduced_y=[];
%  ret.extended_coint_x=[];
%  ret.FF=[];
%  ret.IRF=[];
%  ret.IRF_US=[];
%  ret.IRF_UK=[];
%  ret.D=[];
%  ret.shocks=[];
%  ret.shocks_us=[];
%  ret.p=[];
%  ret.non_nor_shocks=[];


% im=ret.IRF_US(var);
% for i=1:h+1
% hall{i}=repmat(im{1}(i,:),rep,1);
% end
% for i=1:h+1
% hall_{i}=R_{i}(:,1:7)-hall{i};
% end
% for i=1:h+1
% Q{i}=[quantile(hall_{i},[quan,1-quan])];
% end
% for i=1:h+1
% upper_Q{i}=[Q{i}(2,:)];
% end
% upper_bound=reshape(cell2mat(upper_Q'),h+1,7);
% for i=1:h+1
% lower_Q{i}=[Q{i}(1,:)];
% end
% lower_bound=reshape(cell2mat(lower_Q'),h+1,7);
% up_hall=cell2mat(im)-lower_bound;low_hall=cell2mat(im)-upper_bound;
% plot([upper_bound(:,5) ret.IRF_US{var}(:,5) lower_bound(:,5)])
% clear lower_Q upper_Q ret_1 ret_2 ret_3 Q RES R_ res data_ dat u sig_y1 sig_x1 sig* res* E* G*
% upper and lower bounds givethe upper and lower responses of variable var to the other shocks
% ret_4.Q=Q;
% R_=200 cell in which one the irf of variable v to the other 10 shocks in each realizations
% R_iq, i stands for variable, q stands for period
% IT IS POSIIBLE TO USE TO COMMANDS ABOVE TO ATTAIN THE UPPER AND LOWER BOUNDS FOR ANY VARIABLE

% for i=1:rep
% res(:,:,i)=[w([u(:,i)],:)];
% end
% >> for i=1:p
% Y(:,:,i)=[data(1+p-i:end-i,:)];
% end
% >> Y=reshape(Y,T-p,p*k);
% for i=1:rep
% U(:,:,i)=[res(:,:,i) zeros(T-p,(p-1)*k)];
%  end